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Abstract 

We investigate the properties of two standard energy estimators 
used in path-integral Monte Carlo simulations. By disentangling the 
variance of the estimators and their autocorrelation times we anal- 
yse the dependence of the performance on the update algorithm and 
present a detailed comparison of refined update schemes such as multi- 
grid and staging techniques. We show that a proper combination of 
the two estimators leads to a further reduction of the statistical error 
of the estimated energy with respect to the better of the two without 
extra cost. 



1 Introduction 



A detailed understanding of the statistical properties of many-particle quan- 
tum systems is among the most challenging objectives in condensed matter 
physics and physical chemistry. Apart from a few simple model systems, 
analytical approaches can usually only provide an approximative description 
whose accuracy is difficult to control inherently. Computer simulations, on 
the other hand, in principle yield exact results even for complicated systems 
such as, for example, non-linear quantum chains or quantum crystals. A com- 
bination of the two complementary approaches may thus result in valuable 
insights into the physics of quantum systems. 

While in principle an exact method, the difficulty of computer simulations 
is to actually achieve the desired accuracy in practice. For path-integral 
Monte Carlo (PIMC) methods, which we shall consider in this paper, the 
draw-backs are well known. 1 Being a stochastic method, all results are sub- 
ject to statistical errors which, in principle, can be made as small as desired 
by increasing the simulation time. The necessary discretization of the path 
integral, however, requires an extrapolation to the continuum limit where 
standard local update algorithms exhibit a severe slowing down. By this one 
means that successively generated paths, or more generally, configurations 
of a many-particle system, are highly correlated in the Monte Carlo process. 
This effect greatly enhances the statistical errors of the measurements in a 
given computer time. 

More precisely the statistical uncertainty for the mean value O = (1/N m ) 
X^i=i Oi of an observable O = (O) measured in an importance sampling 
Monte Carlo process is given, in general, by the error estimate 



where a 2 . is the variance of the individual autocorrelated measurements Oi 
at "time" i, r int is the integrated autocorrelation time, and N m is the total 
number of measurements. As it appears, there are three ways to reduce the 
statistical error. The most obvious but most expensive one is to increase 
the number of measurements. Since, however, l/y/N m is a rather slowly 
decreasing function and, in the physically interesting continuum limit, it 
may happen that both the variance and the autocorrelation time diverge 
with some power of the discretization parameter, it is much more promising 
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to ask whether the statistical error may also be reduced by constructing 
an estimator for O with a smaller variance or by employing refined update 
algorithms with smaller autocorrelation times. As we shall see below, the 
latter two strategies are to some extent intertwined, a fact which calls for a 
careful analysis. 

In this paper we will focus on energy estimation. It is well known that in 
path-integral Monte Carlo simulations the energy may be measured using two 
different estimators. One is derived by a straightforward differentiation of 
the partition function and will be called for the sake of brevity the "kinetic" 
estimator since it involves an explicit measurement of the kinetic part of the 
energy. 2 The other is based on the virial theorem for path integrals and will 
hence be referred to as the "virial" estimator. 3 ' 4 The "virial" estimator is 
often judged to be the "better" estimator because in the continuum limit its 
variance is much smaller than that of the "kinetic" estimator. 

Early investigations of the "kinetic" and the "virial" estimators focussed 
on their variances. 4 ' 5 In the following years it was pointed out that a correct 
assessment of the accuracy also has to take into account the autocorrela- 
tions, and it was demonstrated that for a standard Metropolis simulation of 
the harmonic oscillator the allegedly less successful "kinetic" estimator gave 
smaller errors than the "virial" estimator. 6 In a further investigation 7 it was 
shown that conclusions about the accuracy also depend on the particular 
Monte Carlo update algorithm at hand since modifications of the update 
scheme such as inclusion of collective moves of the whole path affect the 
autocorrelations of the two estimators in a different way. 

For a fair judgement of the performance of an estimator, one thus has 
to take into account both the variances and the associated autocorrelation 
times. In this paper we analyze this problem for the two standard estima- 
tors in combination with four different update algorithms, namely (1) the 
standard Metropolis algorithm, (2) the multigrid V-cycle, (3) the multigrid 
W-cycle, and (4) the staging algorithm. The Metropolis algorithm is based on 
local moves and exhibits severe slowing down in the continuum limit. 8 ' 9 The 
other three update algorithms employ non-local moves and therefore reduce 
(V-cycle) or even completely overcome (W-cycle, staging) the slowing-down 
problem. 9,10 As our main result we will then demonstrate how the estimation 
of energy can be improved at practically no extra cost by taking a suitable 
linear combination of the two estimators. We shall see that the optimal com- 
bination has to take into account not only the variances of the two individual 
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estimators but also their covariance, i.e., the cross-correlations between the 
two estimators. 

The outline of the paper is as follows. In the next section we first recall 
the definition of the kinetic and virial estimators for the energy and discuss 
some of their basic properties. We then introduce the "combined" estimator 
and present a theoretical analysis of its properties. In Sec. 3 we define au- 
tocorrelation times and describe how our error analysis was performed. The 
various update algorithms are briefly summarized in Sec. 4. The results of 
the numerical simulations are contained in Sec. 5, and in Sec. 6 we close with 
our conclusions and a few final remarks. 



2 Energy estimators 

For didactic reasons we shall illustrate the improved energy estimation for 
simple one-particle quantum systems governed by a Hamiltonian 11 



6 = + ni) - 
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where m is the mass of the particle and V(x) a potential to be specified 
below. Our theoretical considerations are, however, completely general and 
without any additional problems applicable to many-body quantum systems 
as well. When coupled to a heat-bath at inverse temperature (3 = l/Zc^T, 
the canonical partition function is given in the operator representation by 



Z{(3) = Tre-P* = ]T e 



(3) 



where E n are the energy eigenvalues associated with (2). The equivalent 
path-integral representation 12 we wish to simulate reads (K = 1) 



m 



Z((3) = J Vxexp^- jT dr 

where Z^ is the discretized path integral, 

dx k 



jx 2 + V(x(r)) 
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with e = j3/L being the usual discretization parameter. The trace in (3) 
implies periodic boundary conditions, i.e., in (5) we take xo = xl- 

For the potential V(x) we chose two characteristic shapes covering a wide 
range of physical phenomena. The first one is an anharmonic convex poten- 
tial, 

V{x) = 0.5a; 2 + x 4 (CP), (6) 

relevant for studying fluctuations around a unique minimum, and the second 
one is a double-well potential, 

V(x) = -0.5a; 2 + 0.04a; 4 (DW), (7) 

which exhibits tunneling phenomena. 

2.1 The "kinetic" estimator 

From the partition function (3) it is clear that the average energy is defined 
by 

r)\n 7 ~ - V E p-P En 

u = {H} = "W- = M «-*7* = ^0mr- ( 8 ) 

Since d/d(3 = (l/L)d/de it is easy to see that we get, in the path-integral 
representation, U = lini^oo with = (£4) and 

tt L m { x k - x k ^\ 2 li, , 



2/5 2L fc tlV e ) L k=l 

and (...) now means expectation values with respect to the discretized path 
integral (5), which are therefore also L-dependent. The term L/2(3 stems 
from the /^-dependence of the measure. For the sake of brevity we will in 
the sequel refer to the energy estimator (9) as the kinetic estimator 2 of the 
energy identified by a subscript k since it involves an explicit measurement 
of (p 2 /2m) = lim^^Tk)^), with 

being the path-integral estimator of the kinetic part of the energy. Notice 
that the kinetic part is not given by 

m 2 m I ^ (x k - x k ^\ 2 
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as one might have naively guessed, since for all reasonably well behaved 
potentials m{v%)/2 = L/2/3 + 0(1) diverges as L — > oo (and f3 fixed as 
usual). For the proper definition, however, this piece is just needed to cancel 
the divergence coming from the measure. 

Numerically it appears questionable, however, at first sight whether this 
estimator is very useful. The reason is that most of the signal of ((x k — Xk-i) 2 ) 
is used to cancel the explicit L/2(3 term. We would therefore expect that the 
variance of the "kinetic" estimator, 

4 = {ul) - (U k ) 2 , (12) 

diverges linearly with L, and we would hence expect large statistical errors 
in the continuum limit of the path integral even for update algorithms which 
overcome the slowing down problem mentioned in the Introduction. This is 
what we indeed observed empirically; see the discussion in section 5 below. 
What we need is another estimator T for the kinetic energy with the same 
mean, (T) = (p 2 /2m), but smaller variance. In the statistical literature this 
would run under the label of a "variance reduced" estimator. One possible 
candidate is suggested by means of the path-integral version of the virial 
theorem. 



2.2 The "virial" estimator 

In operator language the virial theorem states that 14 
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where we have used the abbreviation V'(x) = dV/dx. From the above dis- 
cussion one may guess that the path-integral analog should read 
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This relation indeed follows from the invariance of the partition function 
under a rescaling of the "field" x k . To show this, we start from eq. (5) and 
rescale the "field" x k — > Xx k , which gives 
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By choosing A 2 = e/m the factor in front of the kinetic term becomes 1/2 and 
also the /3-dependence of the measure vanishes. Differentiation with respect 
to (3 now only affects the potential energy term. Noting that dV(\x)/d\ = 
xdV(y) I dy\ y= xx, and rescaling variables back after differentiation, one readily 
arrives at = (U v ) with the alternative "virial" estimator 3 ' 4 

U, = ^jz XkV\x k ) + -j-]T V{x k ). (16) 

k=l k=l 

From here we deduce the virial theorem (14) with the well-behaved kinetic 
energy estimator T v = ^Y^k=i x k^'( x k)- The subscript v is a reminder of 
the fact that this energy estimator was derived using the virial theorem or, 
alternatively, that, in contrast to the "kinetic" estimator, it only refers to 
the potential V(xk). An alternative derivation of this estimator which makes 
use of the very same rescaling invariance would again start by rescaling the 
"field" Xk in the partition function (5) by x k — > Xxk- Differentiating In Z^ 
with respect to A and putting A = 1, one obtains in fact directly relation 
(14). The virial estimator can also be derived in yet another way. 4 

Notice that using the (incorrect) estimator identifying it (incor- 

-2 

rectly) with (■§—), applying the (correct) identity (operator formulation of the 
virial theorem) , = (^xV'(x)), and identifying (|xV(x)) (correctly) with 
{\%kV (xk)) one arrives (accidentally) at the same result, but the derivation 
is clearly dubious. 

2.3 Discussion of the two estimators 

As explained above, the variance of the measured energies is a property of 
the estimators we employ. Regarding the dependence on the discretization, 
we expect asymptotically for large L a linear divergence of the variance for 
the "kinetic" estimator whereas the variance should stay roughly constant 
for the "virial" estimator. We emphasize that this dependence is completely 
independent of the update algorithm. 

For purposes of illustration we show in Fig. 1 the variance of the two 
estimators as a function of the number of variables L. Since the update 
algorithm only affects the autocorrelation times and is a priori irrelevant 
for the variance of the individual measurements we may choose data for any 
update algorithm. In Fig. 1 we have used the data obtained by the multigrid 
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Figure 1: Variance of the individual energy measurements using the "kinetic" 
and the "virial" estimators for the two potentials (6) and (7) at (3 = 10. While 
the variance of the virial estimator is roughly constant in the continuum limit 
L — > oo, the variance of the kinetic estimator asymptotically diverges as 
a\ = L/2/3 2 . 



W-cycle, cf. the discussion in sections 4.2 and 5.1 and Table 3 below. As 
expected, the measured variances for the virial estimator are indeed roughly 
constant for all values of L. Minor deviations from the asymptotic value of 
al = 0.1348(12) (CP) resp. 0.3564(24) (DW) (cp. Table 3) are observed only 
for the coarsest discretization of L = 8. 

The behaviour is clearly different for the kinetic energy estimator. We 
first observe that for both potentials the divergence turns indeed out to be 
linear for large enough values of L. The straight line in Fig. 1 is the expected 
asymptotic behaviour of a\ — > L/2[3 2 , which in fact fits the data well for 
L = 256, 512, and 1024. In order to account for the discrepancies which are 
observed for smaller values of L, we looked at the first correction term for 
<7k . A straightforward calculation shows that the variance (12) of the kinetic 
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estimator is given quite generally by 

°^w + w\ (x > v ' [Xi)) < (17) 

where C = —/3 2 dU/dj3 is the specific heat. Instead of evaluating C directly 
we first calculate the variance of the virial estimator, 

o* = ^ + ^(lx i V'(z i ) + yV''(x i )) (18) 

which allows us to express a\ as 

4 = ^ + °l - -^V\ Xl ) + \x 2 V"( Xi )). (19) 

For the quartic potential V(x) = uo 2 x 2 /2 + gx 4 we then find 

(IxiVfr) + \x 2 V"{ Xi )) = 2u 2 (x 2 ) + 10g(xt). (20) 

The expectation value on the right hand side can be expanded as a + ai/L 2 + 
a,2 /L 4 + . . .. In order to compute the first correction term ao we now use 
expectation values for the moments for L = 128, assuming that the 1/L 2 
corrections are already negligible for this value of L. At /3 — 10 we obtained 
(in a separate simulation using the staging algorithm) (xf) = 0.25660(46) 
and (xf) = 0.18085(62) for the convex potential (CP, cu 2 — 1, g — 1), resp. 
5.391(18) and 37.33(19) for the double-well potential (DW, u 2 = -1, g = 
0.04). Numerically, we therefore find the correction terms to be ao = 2.3217 
(CP) resp. 4.1500 (DW), and together with the asymptotic values of a 2 = 
0.1348 (CP) resp. 0.3564 (DW) we finally arrive at a 2 - L/2/3 2 = -0.09737 
(CP) resp. —0.0586 (DW). As can be seen from Fig. 1, taking into account 
this first correction does indeed reproduce the data down to at least L = 64. 
Clearly, for smaller values of L we would need to evaluate the higher-order 
coefficients a±, a2, etc., in order to reproduce the values for a\. 

In summary, we have illustrated the expected behaviour that the variances 
of the measurements are roughly independent from L for the virial estimator 
but diverge linearly for large values of L for the kinetic estimator. Before 
going on we finally remark that the two different energy estimators also entail 
two different estimators for the specific heat by virtue of the relation 

<a, v > = -/3 2 ^- (21) 
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2.4 The "combined" estimator 

In the Monte Carlo process we are measuring the energy using either the 
"kinetic" estimator (9) or the "virial" estimator (16) for the configurations 
{xj}i. This procedure yields the two mean values £/ k = J2i=i U^i/N m and 
U y = Yh=i U v ,i/N m as stochastic variables with the same expectation value, 
(C/k) = (U v ) = U, but different statistical errors. Let us here abbreviate the 
variances of the mean values (i.e., the squared errors) by A 2 = (A£/ k ) 2 = 
(Vl) - (Z7 k > 2 and A 2 = (AU V ) 2 = (Vl) - (Z7 V ) 2 . 

It is clear that any linear combination of the form 

U c (a) =aU k + (l- a)U v (22) 

would also give the same expectation value, (U c ) = U, as the individual 
estimators alone. The variance A 2 C = (AU C ) 2 of the combined estimator U c 
would be given by 

A c 2 = ((aZ7 k + (l-a)F v ) 2 )-(aZ7 k + (l-a)I7 v ) 2 

= a 2 A 2 k + 2a(l-a)Al + (l-a) 2 A 2 v , (23) 

where A 2 V = (U\JJ V ) — (U^)(U W ) is the covariance of the mean values of the 
two estimators. Minimizing the variance A 2 with respect to the parameter 
a yields the optimal a as 

A 2 - A 2 

Q ° pt = A 2 + A 2 -2A 2 V - (24) 

Inserting this optimal a opt into eq. (23) one obtains an expression for the 
minimal variance of U c (a) which reads 

A2 = A 2 A 2 - (Ap 2 
copt - a 2 + A 2 -2A 2 v 

1 _ r? 

(25) 



l/A 2 + l/A 2 -2p/A k A v ' 

where p = A 2 v /(A k A v ) is the correlation coefficient of £/ k and U v . If the two 
estimators were completely decorrelated, p = 0, the optimal combination of 
the two estimators would simply be the error weighted mean of the individual 
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estimators. If the variances A£ and A^ furthermore were equal, the error A c 
would be reduced by a factor of \/2, i.e., one would effectively gain a factor 
of 2 in the run time. If on the other hand the two variances were very much 
different we would not gain very much, since the less accurate estimator 
would hardly get any weight in (22) (i.e., a^0ora~l). 

The optimal combination of the two estimators involves, however, the 
covariance of the two estimators and therefore is in general not simply given 
as the error weighted mean of the individual estimators. In fact, erroneously 
ignoring the covariance of the two estimators might lead to the irritating 
result that the variance of the combined estimator is not even reduced with 
respect to the smaller one of the two individual estimators. 

As it turns out, taking properly into account the covariance of the two 
estimators it may even happen that a > 1 or a < 0, and it may also happen 
that the error of the combined estimator is even more reduced than by the 
factor of y/2, which might naively be expected to be an upper bound for the 
gain in accuracy. 

For further discussion, let us assume without loss of generality that A v < 
Ak and let us introduce the ratio k = A v /Ak < 1. The optimal parameter 
a op t may then be expressed as 

k 2 — pn 

^ = 1 + «*-2 P k' (26) 
and the reduction of the error may be judged by looking at the quantity 

r> — ^c.opt _ 1 ~ P /r) 7 \ 

which by definition satisfies < R < 1. The smaller R, the greater the gain 
by using the combined estimator. For an overview of possible situations, we 
plot a op t and R as a function of the ratio n and the correlation coefficient p 
in Figs. 2 and 3. 

In Fig. 2 we first notice that for vanishing covariance, p = 0, the optimal 
interpolation parameter a opt is always in the range < a opt < 1. For com- 
pletely correlated data, p = 1, on the other hand, the optimal interpolation 
parameter is always negative, a opt — k/{k — 1) < 0. 

Taking a look at the error reduction factor R in Fig. 3 we notice again, 
that for completely decorrelated data p — 0, nothing spectacular happens. 
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Figure 2: The optimal interpolation parameter a opt as given in eq. (26) 
as a function of the ratio k = A v /Ak and the correlation coefficient p = 
AL/(A k A v ). 

The gain is best if the variances of the two estimators are equal, k — 1, and 
there is no gain at all, if the smaller variance is negligible compared to the 
larger one, k = 0. The situation is very different though for highly correlated 
data, p — > 1. We first observe that in this case the gain can be much more 
profitable than the best gain of R = 0.5 for completely decorrelated data. 
Assuming, for example, a correlation of p = 0.9 and a ratio of k — 0.5 we 
find that the reduction factor is R w 0.21, i.e. the variance of the combined 
estimator is then roughly 5 times smaller than the smaller of the individual 
variances. It should be noted that for cross-correlated data there is in fact no 
limit to the gain in efficiency, and that, as a general rule, the gain is largest for 
two estimators which are highly correlated and have very different variances. 
Notice that for p ss the inverse gain factor R is a decreasing function of k 
(highest gain for comparable variances), while for p 1 the situation is just 
reversed, i.e., R increases with k (highest gain for very different variances). 
In order to illustrate the combination of the two estimators with actual 
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Figure 3: The reduction factor R as given in eq. (27) as a function of the 
ratio k = A v /Ak and the correlation coefficient p = A^ v /(AkA v ). 

simulation data (using data for the convex potential obtained by multigrid W- 
cycle updating), we show in Fig. 4 the relative statistical error as a function 
of the interpolation parameter a. In order to simulate the measurement of an 
arbitrary combination of the two estimators, the data points in Fig. 4 were 
obtained from a times series of an estimator U c (a) for arbitrary a which was 
generated from the time series of each of the individual measurements of 
and U v . The time series for U c (a) was then subsequently analysed for each 
a by the usual jackknife blocking procedure to obtain the error estimate (cf. 
section 3.2). The solid lines interpolating these data points, on the other 
hand, were computed using eq. (23) with the variances and covariance of 
and U v obtained by the jackknife method on the basis of 100 blocks. 

We see that the theoretically expected error (23) reproduces indeed the 
empirically measured errors for a linear combination of the two estimators at 
each measurement. We emphasize that the combination can hence always be 
done post simulation without any restriction, as long as the variances and the 
covariance of the two estimators have been measured (i.e., it is not necessary 
to store the complete time series of and U v which, in some applications, 
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Figure 4: Relative statistical error AU C /U C of the linear combination U c of 
the two energy estimators for the convex potential as a function of the inter- 
polation parameter a (using the W-cycle update algorithm). Data symbols 
denote jackknife averages over 100 blocks, and the solid lines are computed 
according to the theoretical prediction (23), using as input only the (jack- 
knife) variances of the virial and the kinetic estimator (corresponding to 
a — 0, resp. a — 1) and the (jackknife) covariance of the two estimators. 
The arrows indicate the values of optimal a according to eq. (24). 

might cause disk-space problems). 

The arrows indicate the locations of the optimal values of a computed 
according to eq. (24). These values are also reported in Table 1 together 
with the correlation coefficient p and the measured energies using the optimal 
combination. We note that for the coarsest discretization of L = 8 we do 
indeed find a case where the optimal interpolation parameter a falls outside 
the usually expected range of [0,1]. Unfortunately, in this application the 
theoretically possible error reduction for correlated estimators is not realized, 
and the physically interesting data are those for large L, where the correlation 
coefficient is only about p = 0.3 — 0.5. Very similar plots were obtained for 
the case of the double-well potential. 
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At this point it should be stressed that, in general, the question of how 
much the error can be reduced in the simulation of a particular problem at 
hand appears to be an empirical question. As a general lesson, we emphasize 
that the fact that two independent estimators are correlated does not neces- 
sarily mean bad news. The gain in efficiency crucially depends on both the 
correlation and the ratio of the individual variances. A further discussion 
will be given below in section 5.4. 

3 Errors in the Monte Carlo process 

In this paper we investigate both the variances and the autocorrelation times 
associated with different energy estimators. While the variance is a property 
of the partition function and the estimator alone, the autocorrelation time 
also depends on the update scheme and represents the relevant quantitative 
measure for the dynamics of the Monte Carlo process. 

In general, the autocorrelation time r is proportional to some power of 
the correlation length £ of the system, r oc f , with a dynamical critical 
exponent z. For standard local update algorithms the exponent z is close to 
2 as can be argued in a simple random walk picture. For spin systems or field 
theoretic models undergoing a second-order phase transition the correlation 
length diverges and as soon as £ exceeds the linear size L of the lattice the 
autocorrelation times diverge like r oc L z . In statistical mechanics and lattice 
field theory this problem is known under the name of critical slowing down. 

In path-integral simulations a very similar problem occurs in the contin- 
uum limit e — > with (3 fixed, or, equivalently, L — > oo. The reason for this 
slowing down is easily understood. The correlations (xkXk+i) only depend on 
(3 and on the gaps between the energy levels. Hence the correlation length £ 
only depends on the physical parameters at hand, and consequently always 
diverges linearly with L if measured in units of the lattice spacing e. Thus 
we expect that the autocorrelation time for path-integral simulations grows 
as 

r oc L\ (28) 

with a dynamical critical exponent z as well, and that for standard local 
algorithms z — 2. Note that in contrast to the infinite- volume limit in sta- 
tistical mechanics, for path integrals also the Hamiltonian, i.e., the exponent 
in (5) changes its form in the limit L — > oo which is a characteristic feature 
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of path-integral simulations. As a consequence slowing down occurs in the 
continuum limit for any fixed (3 and any set of potential parameters. 



3.1 Definition and measurement of autocorrelation times 

Since the measurement of autocorrelation times will be of central concern in 
the remainder of this paper let us briefly recall the basic definitions. 17,18 

In general, if Oi denotes the ith. measurement of an observable O in 
the Monte Carlo process the (normalized) autocorrelation function A(j) is 
defined by 

_ (O t O l+J ) - (Of 

m ~ (of)-(o^ ■ (29) 

In Monte Carlo simulations we are always dealing with a finite number of 
measurements N m . Nevertheless, computing the variance for the mean O 
of measurements is straightforward and yields an error estimate of the form 
AO = \j2Ti nt (N m ) ^Ja 2 /N m , where N m is the number of measurements used 
to compute the mean value O. Here T int (k) is an integrated autocorrelation 
time given by 19 



1 



k 



2 i=l 



= £ + E^0') 1-tt • (30) 



N„ 



The effective statistics is thus reduced to N e g = N m /2T m t(N m ). Or, in other 
words, to achieve a given error AO the run-time (i.e. the budget) has to be 
increased by a factor of 2r int (iV m ,). 

For large j the autocorrelation function A(j) usually decays like an ex- 
ponential 

A(j) ^ a exp(-j/r cxp ), (31) 

where r cxp denotes the exponential autocorrelation time and a is some con- 
stant. Since, in general, all these quantities depend on the observable under 
consideration we will indicate the relevant observable by an additional sub- 
script unless it is clear from the context which observable is meant. 

Due to the usual exponential decay of the autocorrelation function A(j) 
and the large number of measurements N m in a Monte Carlo simulation, the 
factor 1 — ji— in (30) accounting for the finite size of the statistical sample 
may in practical applications savely be neglected. In practice, one therefore 
usually computes the autocorrelation time by computing the sum in eq. (30) 
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without the correction factor and by cutting it off self-consistently at some 
k max such that n cut T int (k max ) « k max « N m . Usually n cut is set equal to 
^cut = 6 or 8. As long as integrated and exponential autocorrelation times are 
roughly the same this method gives reliable estimates for r int . If exponential 
and integrated autocorrelation times are very much different from each other 
this method tends to underestimate the integrated autocorrelation time. This 
typically happens when the faster decaying modes neglected in (31) are still 
important for relatively large time lags j. 

A more reliable way of determining the autocorrelation times in this case 
is to rewrite the integrated autocorrelation time as 20 

oo 

Tint(fc) = r int -a ^ exp ( - j / T cxp ) (32) 
j=k+l 

= Ti nt -a 6XP / ^ -, exp{-A;/r eX p}, (33) 
1 - exp{-l/r cxp } 

where Tint = T int (oo). A three-parameter fit of r int (A;) then gives both the 
integrated and the exponential autocorrelation times Tint an d T cxp simultane- 
ously which is a further advantage of this method. 

In Fig. 5 the behaviour of A(j) and r- int {k) of the virial estimator is il- 
lustrated for one typical example (CP potential, V-cycle multigrid update, 
(3 — 10, L — 512). A two-parameter fit of the autocorrelation function A(j) 
according to eq. (31) in the range j = 5, ... ,15 for this case gave values of 
a = 0.619(33) and r cxp v = 7.23(40). A three-parameter fit of the integrated 
autocorrelation time according to eq. (33) in the same range yielded values 
of a = 0.614(42), r cxPjV = 7.26(48), and T intjV = 4.98(15). A determination 
of the integrated autocorrelation time via a self-consistent cutoff at 8Ti n t, v on 
the other hand yielded a value of r int v = 5.07(21) (cp. Table 3). Clearly all 
values are consistent within error bars. 

3.2 Blocking and jackknife procedures 

Another way of estimating the true error of the Monte Carlo simulation is 
to divide the N m measurements into riu blocks of size iV bl = N m /nu and 
compute the averages for each block separately. The block averages are then 
again stochastic variables with the same mean but reduced autocorrelations. 
In fact, if the block length is appreciably larger than the integrated auto- 
correlation time for the observable under consideration (20T int , say, would 
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clearly be sufficient) than the block averages are (almost) uncorrelated and 
we can estimate the error of the observable from the variance of the block av- 
erages, AO = y / o"o i bi/ n bi- Since nu should be at least around 20 to allow for 
a statistically meaningful estimation of (Tq U , the number of measurements 
per block is always much smaller than the total number of measurements, 
iVbi <C N m . For observables that can be estimated only with so-called biased 
estimators this may result in severe systematic errors when Nu becomes too 
small. 

It is therefore gratifying that bias-reduced and thus more accurate es- 
timates of the errors can be achieved by using so-called jackknife blocking 
techniques. 21 The main difference is that here the block averages are taken 
over the whole run with only one block excluded. These jackknife block 
averages have consequently a much reduced variance 0"e>j_bi an d the error 
estimate now reads AO = \JJjib\ — l)aQ^_ hl /nu, where the factor (n w — 1) 
corrects for the trivial correlations between different jackknife blocks (because 
each measurement enters in all but one jackknife block - this has nothing to 
do with autocorrelations). It should be stressed that for unbiased estimators, 
such as arithmetic mean values, the blocking and jackknife blocking method 
give identical results. The advantages of jackknife blocking only show up 
for biased estimators, such as those for the specific heat or autocorrelation 
functions. 

Using the blocking or jackknife blocking error estimate the integrated au- 
tocorrelation time for an observable O can then also be obtained by inverting 
the standard error formula as 

(AOfN m 

where AO is the error measured in the blocking analysis and Oq. is the vari- 
ance of the single (unblocked) measurements. It should be stressed that this 
method of estimating the integrated autocorrelation time is only valid if the 
block length iVbi is several times larger than the true integrated autocorre- 
lation time. Otherwise T mt will be systematically underestimated. At the 
same time we want to have, of course, as many blocks as possible in order 
to reduce the statistical error on the error estimates. Since ribiiVbi = N m is 
fixed for a given statistics, these are conflicting requirements. If the mea- 
surements are not too strongly correlated, a reasonable compromise is the 
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choice ribi ~ Nu ~ y/N, 



1.000 



0.100 



AG) 



0.010 



x int (k) 3 




CP, V-cycle 
virial estimator 

L = 512 



(a) 



10 20 
j 



30 




CP, V-cycle 
virial estimator 

L = 512 



(b) 



20 



40 



60 



Figure 5: (a) The autocorrelation function A(j) on a logarithmic scale and (b) 
the integrated autocorrelation time Ti nt (k) of the virial estimator as obtained 
with the V-cycle multigrid algorithm for the convex potential (CP) at /3 — 10 
and L = 512. Solid lines show fits according to eq. (31) resp. (33). The 
asymptotic value of T int (A:) quoted in Table 3 is r int)V = 5.07(21). 
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4 Update algorithms 



Before presenting our results in the next section, we will briefly review the 
update algorithms we have investigated, notably the multigrid method and 
the staging algorithm. 

4.1 Local algorithms 

The most extensively used algorithm for PIMC simulations still is the stan- 
dard Metropolis algorithm. 22 As far as accuracy and efficiency are concerned 
it has the well-known serious drawback of local algorithms that in the con- 
tinuum limit e — j3/L — > the autocorrelation times diverge quadratically 
with the grid size, r oc L 2 , resulting in a severe slowing down of the Monte 
Carlo process, see the discussion in section 3. In spite of this drawback it is 
still widely used for its simplicity even though a lot of computer time might 
be saved by using more refined algorithms. 

A local algorithm which does reduce slowing down to some extent is the 
hybrid overrelaxation method. 23 This algorithm mixes (deterministic) over- 
relaxed updates of the path with stochastic Metropolis sweeps. If the ratio 
of overrelaxed and Metropolis sweeps is chosen properly, i.e., in proportion 
to the spatial correlation length, this algorithm reduces slowing down to a 
linear divergence, 8 ' 24 r oc L. 

More successful are non-local algorithms which will be briefly described 
in the next two subsections. 

4.2 Multigrid method 

The basic idea of the multigrid approach 25 ' 26 is to perform non-local updates 
of the variables by working on a set of successively coarser discretizations of 
the time axis ("grids") in order to take into account long wavelength fluc- 
tuations of the paths more efficiently. To this end one performs collective 
updates on different length scales by visiting various coarsened grids in a 
systematic order as extensively discussed in the context of partial differential 
equations. 26 The auxiliary variables on the coarsened grids are then interpo- 
lated back to the finer grids and eventually to the original grid using some 
specific interpolation scheme in a recursive manner. Equivalently, one may 
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also view the multigrid approach from a unigrid point of view where the up- 
date on a coarsened grid corresponds to a simultaneous move of a group of 
adjacent variables on the original grid. Using the so-called piecewise constant 
interpolation scheme for example, this amounts, in the unigrid viewpoint, to 
proposing moves for blocks of 1, 2 d , 4 d , . . . , V = L d = 2 nd adjacent variables 
in conjunction. Here d is the dimension of the system under consideration. 
For PIMC simulations we usually have d — 1, even for quantum chains 
or quantum crystals since it may well suffice to use multigrid acceleration 
only along the Trotter direction of the discretized path integral. 27 Partic- 
ular successful sequences of length scales 2 k are the so-called V-cycle with 
k = 0, 1, . . . , n — 1, n, n — 1, . . . , 1, 0, and the W-cycle whose graphical rep- 
resentation looks like the letter W (for n = 3, e.g., this is k = 0, 1, 2, 3, 2, 
3, 2, 1, 2, 3, 2, 3, 2, 1, 0). 28 The update at level k thus consists in consider- 
ing a common move Ax for all 2 kd variables of one block, x^ — > x^ + Ax, 
i G block, computing the associated energy change and applying the usual 
accept /reject criterion. 

The multigrid approach thus has a number of parameters which may be 
adjusted to suit the problem at hand, notably the choice of the interpola- 
tion scheme, the block length, the recursion scheme (e.g. V- or W-cycle), 
the number of updates on each grid in going down the recursion scheme 
(presweeps) and in going up again (postsweeps), and the Metropolis param- 
eters for the coarsened grid updates. In the Gaussian case it is known that 
by using the piecewise constant interpolation scheme the V-cycle leads to a 
linear divergence of autocorrelation times, r oc L, while the more successful 
W-cycle beats slowing down completely, r oc const. 9 

Let us finally emphasize that the multigrid method can easily be applied 
to general d- dimensional lattice field theories. In fact, it is in this context 
where the stochastic multigrid Monte Carlo formulation appeared first in 
the literature, 25 and only quite recently these ideas have been adapted to 
path-integral simulations. 8,9 ' 29 

4.3 Staging algorithm 

The basic idea of the staging algorithm 30 is to rewrite the discretized quan- 
tum statistical partition function (5) in such a way that a sequence of j 
adjacent variables can be updated one by one independently, i.e., by effec- 
tively removing the coupling between the variables. The coupling between 
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variables stems from the kinetic term and, in the continuum limit, e — > 0, 
the kinetic term mj (2e)(xk — Xk-i) 2 numerically starts to dominate over the 
potential energy term eV(x k ). The kinetic term for adjacent variables can, 
however, trivially be rewritten as 



(X k +1 - X k ) 2 + (x k+2 - X k+ if = (Xfc+2 - x k ) 2 /2 + 2(x k+1 



X 



\ 2 

k+l) i 



(35) 



with xl +1 = (xk + x k+2 )/2. The variable x k+ \ can thus be decoupled from 
its neighbours. The crucial observation for the staging algorithm is that this 
can be iterated as 



3j~j~y ( x k+j+i~ x k) 2 + ^— {x k+ j — ) 2 > 

(36) 

x k + Xk+j+i)/(j + 1)- A sequence of free particle 
propagators can thus be rewritten as 30 



— (Xk+j— x k) 2 ~\~ { x k+j+l ~ x k+j) 2 



with new variables x* k+ j 
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where renormalized masses m 8 = -^m have been introduced. Selecting the 
end points x k and x k +j+i of some segment of the discretized path with j 
"beads" in between, one can thus perform a (recursive) change of variables 
x k+ i — > x* k+i) i = 1, . . .j, in the discretized partition function (5). For the 
staging segment, this would eliminate the nearest neighbour coupling stem- 
ming from the kinetic energy. For the variables of the staging segment the 
partition function hence reduces to a collection of independent oscillators 
moving in an external potential which depends on the transformation of the 
variables. The staging variables may then be updated using Gaussian random 



variables u k = x k 



with "masses" m^, and a Metropolis like acceptance 



rule for the external potential. 

In contrast to the multigrid method the staging algorithm only allows for 
one single tunable parameter, the length j of the staging segment. The opti- 
mal choice of the staging parameter j depends on the observable of interest, 
its estimator and on the discretization. But since it is the greater mobility 
of the variables of the staging segment which reduces autocorrelation times 
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the optimal staging length j opt scales with the correlation length along the 
path, i.e., it scales with the number of variables L in the discretized partition 
function. 10 For an appropriate choice of j opt the staging algorithm then com- 
pletely reduces slowing down in the continuum limit, r oc const. Thus, both 
the staging and multigrid W-cycle represent refined PIMC update schemes 
which beat the continuum slowing down. Elsewhere we have compared the 
two algorithm and discussed their mutual merit in more detail. 10 

5 Results 

5.1 Simulation details 

For each of the four update algorithms described above (Metropolis, multigrid 
V- and W-cycle, staging algorithm) we simulated the path integral (5) for 
grids of size L = 2 3 = 8 up to L = 2 10 = 1024 sites. In all our simulations 
the mass m was set equal to 1 and the inverse temperature was equal to 
(3 = 10. In our simulations we performed N m = 100 000 updates of the path 
after discarding 5 000 initial updates of the path for thermalization. 

In the case of the Metropolis algorithm an "update of the path" here 
means n e sweeps over the full path with single-hit updates of each site with 
roughly 50% acceptance probability. In order to accurately assess autocor- 
relation times, the parameter n e was adjusted in such a way that the au- 
tocorrelation times in units of measurements were comparable for all grid 
sizes L. For the convex potential this could always be achieved by setting 
n e — 1, except for L = 256, 512, and 1024 where we had n e = 5, 20, and 
80, respectively. For the double- well potential we started out with n e — 1, 
1, 1, and 2 for L = 8, 16, 32, and 64. For the larger grid sizes, however, the 
autocorrelation times in units of single sweeps turned out to be so different 
for the two estimators that we actually performed two sets of simulations 
with different choices of n e (the combined estimator then obviously makes 
no sense). In the first set, focusing on the autocorrelation time for the kinetic 
estimator, we could still use n e — 1, except for the largest grid L = 1024 
where we set n e = 6. In these runs the autocorrelation times for the virial 
estimator turned out to be far too big to be measurable on the larger grids. 
To satisfy our own curiosity and to measure autocorrelation times for the 
virial estimator as well, we therefore performed a second set of simulations 
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with n e = 15, 60, 150, and 600 for L = 128, 256, 512, and 1024, respectively, 
even though this clearly requires employing CPU resources out of proportion. 
In the following, all autocorrelation times for the Metropolis algorithm will 
be given in units of single sweeps. 

In the case of the multigrid algorithm an "update of the path" means a 
complete V- resp. W-cycle with n prc = 1 presweeps and n post = postsweeps. 
On each grid we performed single-hit Metropolis updates with an acceptance 
rate of 40% — 60% which for the system at hand could be achieved with the 
same maximal step width on all grids. 

In the case of the staging algorithm, an "update of the path" means 
int(L/(j — 1)) calls to the staging routine which moves j — 1 adjacent vari- 
ables at each call. The choice of the staging length j is shown in Table 5 
below. For the analysis of the combined estimator we used the same j as for 
the virial estimator. This will be discussed in more detail below in section 
5.4. Notice that the above definition in general implies updates of less than L 
variables. We have therefore rescaled the actually measured autocorrelation 
times by a factor (int(L/(j — l)))/(L/(j — 1)). This enables a direct com- 
parison of the staging algorithm with the standard Metropolis algorithm. A 
comparison with the multigrid method has to take into account a constant 
factor for the V-cycle (which of course depends on the implementation but 
should theoretically be ~ 2) and an extra factor of log L for the W-cycle. 

After each update of the path we measured the internal energy using both 
the "kinetic" estimator (9) and the "virial" estimator (16). The time series 
of these measurements were then analysed by jackkniving the data on the 
basis of 100 blocks. Autocorrelation times were obtained by cutting r int (/c) 
self-consistently at k max = 8. The reported errors for the autocorrelation 
times were again obtained by jackkniving. 

With these remarks in mind we emphasize that our data now allow for 
a precise and detailed comparison of the commonly used energy estimators 
taking into account the dynamics of different update schemes. 

5.2 Results for the energy 

Tables 1 and 2 show the measured energies for the convex potential (CP) 
and for the double-well potential (DW) using the "kinetic" estimator and 
the "virial" estimator as well as using the optimally combined estimator. 
For ease of comparison we also give for each energy measurement the rela- 
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tive error in percent. In addition, for the combined estimator, we list the 
correlation coefficient p and the parameter a opt used for computing the com- 
bined estimator according to eq. (22). 

Since in our simulations we saved the time series of the U^- and U v - 
measurements, we were able to compute a posteriori time series of U c for 
any a which could then be analyzed in the same fashion as the (run-)time 
series for the "kinetic" and the "virial" estimator. The energy values and the 
jackknife errors reported in Tables 1 and 2 were thus obtained by analyzing 
the time series of the combined estimator U c (a opt ) for the optimal choice of 
a. From the discussion in section 2.4, in particular Fig. 4, it should be clear 
that we could equally well have applied error propagation using the variances 
and the covariance of the two original estimators. 

The analysis of a times series of U c has the advantage that, for a more 
detailed discussion, the variance of the individual measurement and the au- 
tocorrelation time of the combined estimator can be computed as well. In 
Tables 3 and 4 we list for each estimator the variance cr^ vc of the individual 
measurements as well as their integrated autocorrelation times Ti nti k,v,c- As 
already mentioned, the reported errors for these quantities were obtained 
using jackknife blocking. 

Looking at the energy values in Tables 1 and 2, the first thing we notice is 
that all three estimators give indeed compatible energy estimates within error 
bounds and all estimates converge to the correct continuum energy of U = 
0.80377 (CP) resp. -0.903965 (DW). These values can be easily obtained by 
numerical integration of the associated Schrodinger equation. For the convex 
potential this is basically the ground-state energy Eq, since already the first 
excited state with E\ = 2.736 is strongly suppressed at (3 = 10 and does not 
affect the significant digits of U . For the double-well potential the value of 
U was obtained by using Eq = -0.913371, E 1 = -0.892348, E 2 = 0.029846, 
and E 3 = 0.37813. 

One also recognizes a distinct finite-size dependence which on theoretical 
grounds should be proportional to L~ 2 . For the most accurate measurements 
in our simulations (using the combined estimator and the W-cycle) we find 
that only for L > 256 (CP) resp. L > 128 (DW) the discretization error is 
no longer relevant. We mention in passing that this systematic error may be 
reduced using the Takahashi-Imada scheme 31 based on higher-order Trotter 
decomposition. Since this modification only amounts to adding a local term 
to the potential energy this scheme is perfectly compatible with all update 
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algorithms discussed here, and as far as the qualitative behaviour of statistical 
errors is concerned we do not expect any significant deviations from the 
conclusions drawn in this study. 

5.3 Autocorrelation times and statistical efficiency 

The integrated autocorrelation times for the two potentials and the four 
update algorithms collected in Tables 3 and 4 are plotted in Figs. 6(a)-(h) 
against L on a double-logarithmic scale. As discussed in section 2.3 and il- 
lustrated in Fig. 1, the variance does not depend on the update algorithm 
but it strongly depends on the estimator. As we will see, the different de- 
pendencies of the variance and the autocorrelation times on the number of 
variables L combine in a peculiar way in the errors A = AU which are plot- 
ted in Figs. 7(a)- (h). Comparing the "kinetic" and the "virial" estimators, 
the combination often leads to a crossover for the final error. The error of 
the combined estimator is, however, always smaller than the better of the 
two estimators. 

Metropolis algorithm: In Figs. 6(a) and (b) we observe for the autocor- 
relation time of the virial estimator the expected quadratic divergence. By 
fitting the power-law ansatz r int = aL z to the data for the convex potential 
(CP) with L = 128, 256, 512, and 1024 we obtain indeed an exponent of 
z = 2.017(13). The autocorrelation times of the kinetic estimator in this 
case behave rather differently, and at first glance it may seem that they do 
not show the expected L 2 -divergence. It is clear, however, that the plotted 
fit with an "effective" exponent of z — 1.372(16) is not really justified since 
there is still distinct curvature for the last three data points (L = 256, 512, 
1024). For the double- well potential the difference between the two esti- 
mators is even more pronounced. For the virial estimator we obtain in the 
range 64 < L < 1024 again an excellent fit with an exponent of z — 1.986(54), 
which is fully consistent with the expected value of two. For the data of the 
kinetic estimator in the same range the quality of the fit is greatly reduced 
and the "effective" exponent of z — 0.461(15) is again found to be much 
smaller. By successively discarding small values of L in the fit we observed a 
definite trend to larger values of z, but the available sizes of L are still much 
too small to see the truly asymptotic behaviour for the kinetic estimator. 
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We may, however, compare the autocorrelation times with those of the 
moments m 2 = Y^f=i x 1 an d = X) i=1 x t reported in an earlier investi- 
gation, 9 which are of the same order of magnitude. The same is true for 
the autocorrelation times of the correlations x^Xk-i which we have elsewhere 
found 24 for the convex potential to increase as 1.078(14), 1.903(22), 4.466(94), 
13.53(70), 52.3(6.6), and 227(21) for lattices of size L = 8, 16, 32, 64, 128, 
and 256. Surprisingly enough, the autocorrelation times of the kinetic esti- 
mator are much smaller and start to diverge only much later. Recall that due 
to the periodic boundary conditions and translational invariance the kinetic 
estimator is a linear combination of TO2, wi4, and xj-Xk-i- The particular 
combination, however, has a much smaller autocorrelation time. 

As a consequence, contrary to previous expectations, for the Metropolis 
algorithm the kinetic estimator turns out to be much more accurate than the 
virial estimator for the parameters under consideration. This is clearly seen 
in Figs. 7(a) and (b). Note, however, that this result depends crucially on the 
observed autocorrelation times in the range of L-values investigated here. If 
also for the kinetic estimator the autocorrelation times asymptotically diverge 
as r intj k oc L 2 , then we should find for very large L that (A£/ k ) 2 oc o-j^mt.k oc 
L 3 , while the error for the virial estimator exhibits a weaker L-dependence 
of (AZ7 V ) 2 oc L 2 . 

Since modified update algorithms such as the staging algorithm and the 
multigrid W-cycle eliminate slowing down the virial estimator for these up- 
date schemes will asymptotically always be more favorable. 

Multigrid V-cycle: The different behaviour of the autocorrelation times 
for the kinetic and virial estimators is also found, albeit less pronounced, 
for the V-cycle. Fitting our data for L = 256, 512, and 1024 in Figs. 6(c) 
and (d) we find for the virial estimator exponents of z — 0.959(54) (CP) 
and z = 0.808(42) (DW), supporting the expected value of z — 1. The 
corresponding exponents for the moments 1712 of z = 0.8356(92) (CP) and 
z = 0.715(27) (DW), obtained elsewhere 9 fitting data for L = 128, 256, 
and 512, are compatible with these values taking into account the fact that 
those data still showed some upward curvature and were hence conjectured 
to underestimate the asymptotic behaviour. 

Again, however, the kinetic estimator has much smaller autocorrelation 
times. For L < 1024 the fits give here effective exponents of z — 0.135(18) 
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(CP) and z = 0.078(23) (DW). Even though we can of course not exclude the 
possibility that the data again start to diverge much faster for larger values 
of L, our data do not show any tendency of upward curvature. 

Taking into account the fact that also the variance for small values of L 
is much smaller for the kinetic estimator we find that, for small L, the actual 
errors for the kinetic estimator are much smaller as well. In Figs. 7(c) and 
(d) we see that with increasing L the difference decreases, and for L m 50 
both estimators would roughly yield the same error. For large values of L the 
virial estimator gets still better but for values in the range 256 < L < 1024 
the errors for the two estimators seem to increase with roughly the same 
effective exponent. 

Multigrid W-cycle: The W-cycle exhibits a considerably improved dy- 
namical behaviour and the magnitude of the autocorrelation times shown 
in Figs. 6(e) and (f) is now greatly reduced to about unity. Notice that in 
contrast to all other update algorithms here the virial estimator has smaller 
autocorrelations than the kinetic estimator. In particular the divergence 
with increasing L is much weaker for both estimators. Only for the convex 
potential we still find a slight increase of the autocorrelation times for the 
virial estimator with an exponent of z — 0.1087(65), while the corresponding 
exponent for the kinetic estimator, z = 0.052(11), is compatible with zero. 
For the double-well potential both exponents [z = 0.040(63) for the virial 
and z = 0.028(12) for the kinetic estimator) are consistent with zero. Hence 
the continuum slowing down problem is solved for the W-cycle. Again these 
estimates are in good agreement with the corresponding exponents for the 
moments m 2 of z = 0.1043(29) (CP) and z = -0.015(11) (DW) obtained in 
our earlier investigation. 9 

Since the autocorrelation times no longer diverge and since the behaviour 
of the variances does not depend on the update algorithm, we expect a 
crossover for the actual errors associated with the virial and the kinetic es- 
timator. Such a crossover was already observed for the V-cycle in Figs. 7(c) 
and (d), but for the W-cycle in Figs. 7(e) and (f) it is much more pronounced. 
Basically this simply reflects the behaviour of the variances a% and a\ shown 
in Fig. 1, which exhibit a clear crossover around L = 64 for the convex and 
around L = 100 for the double-well potential, respectively. Since for large 
values of L the statistical errors of the virial estimator remain roughly con- 
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stant, it here always outperforms the kinetic estimator whose errors increase 
due to the linear L-dependence of the variance. 

Staging algorithm: From Figs. 6(g) and (h) it is obvious that also the 
staging algorithm eliminates slowing down. Notice that for the staging algo- 
rithm the virial estimator has larger autocorrelations than the kinetic esti- 
mator, while for the W-cycle it is just the other way around. 

It should be kept in mind, however, that the staging algorithm requires 
the choice of an optimal staging length j opt . This problem was discussed 
elsewhere 10 for the virial estimator by explicitly looking at the autocorrela- 
tion times as a function of the staging length. It was shown that the optimal 
choice scales with the number of variables L, and that the corresponding op- 
timal acceptance rates stay roughly constant for different L. Here we extend 
this discussion and list in Table 5 the values of j op t and the corresponding 
acceptance rates for both the virial and the kinetic estimator. We see that 
the rule of thumb of some fixed acceptance probability regardless of the esti- 
mators is rather misleading. In fact, while for the virial estimator the optimal 
value of j corresponds to an acceptance rate of about 55%, for the kinetic 
estimator the optimum is at about 85% — 90%. 

It should be observed that for an optimal performance we had to use a 
different j opt for the two estimators, i.e., used data from different simula- 
tions. For the convex potential and the largest grid of L = 1024 the kinetic 
estimator had the lowest autocorrelation time of T°J^ k = 1.097(21) (cp. Table 
3) for a staging length of j op t = 56 with an acceptance rate of 90%. Mea- 
suring the energy, on the other hand, for j = 176 with an acceptance rate of 
56%, as was best for the virial estimator, produced an autocorrelation time 
of r mt,k — 1.483(37). The virial estimator, on the other hand, had its lowest 
autocorrelation time of t°J^ v = 2.48(56) for j opt = 176. Here the autocor- 
relation time is in fact more than doubled to a value of r inti v = 4.53(17) if 
the optimal staging length j = 56 for the kinetic estimator was used. This 
is nothing but yet another example of the subtle interplay between update 
algorithm and energy estimator which calls for some care when one aims at 
optimizing the efficiency of PIMC simulations. 

A bad choice of the staging length - if only properly scaled with L - only 
affects the absolute value of the autocorrelation times. 10 If the staging length 
increases linearly with L the corresponding autocorrelation times are roughly 
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constant and do not depend on L. This is confirmed by the data displayed 
in Figs. 6(g) and (h) where we plotted fits with exponents of z — 0.008(17) 
(CP) and z = —0.005(20) (DW) for the virial estimator, and with exponents 
of z = -0.012(13) (CP) and z = 0.005(12) (DW) for the kinetic estimator. 
All these exponents are fully consistent with zero and we conclude that the 
staging algorithm eliminates slowing down just as well as the multigrid W- 
cycle. Looking at the error of the energy estimation in Figs. 7(g) and (h), 
which contains also the L-dependence of the variance, we hence observe the 
very same crossover as we did for the W-cycle. 

5.4 The combined estimator 

If the choice between the two estimators were exclusive, one would have to 
conclude that the best estimator would depend on details of the simulation. 
As discussed in section 2.4, however, one can always combine the two esti- 
mators and can hence always obtain errors which are even smaller than the 
value for the better of the two estimators. Looking at the errors displayed in 
Figs. 7(a)- (h), one sees that for the convex potential the gain by combining 
the two estimators is best for very small values of L where a is outside the 
range [0, 1] (cp. Tables 1 and 2). Another situation where the combination 
appreciably reduces the error occurs, for both potentials, if the actual errors 
of the two estimators are of the same magnitude. 

As was already suggested by Fig. 4, for refined update algorithms and 
large values of L the optimal interpolation parameter is always close to 
which reflects the fact that the virial estimator here always wins the race. 
For the Metropolis algorithm where no crossover occurs, the opposite is true. 
Here the optimal a for large L is close to 1. For medium values of L, in 
the crossover region for the V- and W-cycle and the staging algorithm, any 
values of a were obtained. 

Let us take a look at the autocorrelation times for the combined esti- 
mator displayed in Figs. 6(a)-(h). For the Metropolis algorithm the auto- 
correlation times (scaling effectively with exponents z = 1.413(21) (CP) and 
z = 0.469(51) (DW)) are close to but not necessarily smaller than those of the 
kinetic estimator. For the V-cycle the autocorrelation times of the combined 
estimator (z = 0.523(38) (CP) and z = 0.388(33) (DW)) are appreciably 
larger than those of the kinetic estimator but smaller than those of the virial 
estimator. For the W-cycle (z = 0.136(12) (CP) and z = 0.0110(88) (DW)) 
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and for the staging algorithm (z = 0.021(17) (CP) and z = 0.042(17) (DW)) 
the situation is just reversed. Here the combined autocorrelation times be- 
have qualitatively as the ones for the virial estimator. 

Regardless whether the autocorrelation times of the combined estimator 
are close to the smaller ones of either of the original estimators it always 
yields errors which are slightly smaller. The gain is most pronounced for 
very small lattices and in the crossover region. 

The best gain in the reduction of the error by combining the two estima- 
tors was obtained for L = 8 and the Metropolis algorithm where the error 
of the better estimator was in fact reduced by a factor of 2. For L = 64, 
the double well and the V-cycle the errors of 0.31 for the kinetic and 0.30 
for the virial estimator were still reduced by a factor of 1.5 to an error of 
0.21. Note that this reduction would be equivalent to a factor of more than 2 
for the actual computer time needed to obtain the same reduction by better 
statistics for the individual estimators. 

Unfortunately, for most cases the gain by combining the estimators is 
much smaller. Realistically one may expect a gain of, say, 10% in the final 
error by the combination which may seem moderate a gain. It should be 
kept in mind, however, that such a 10% gain in the error corresponds to 
something like a 20% reduction of the computer run time needed to achieve 
the same error by a simply increasing the statistics. And we emphasize again 
that the reduction gained by combination of the estimators may be obtained 
after the simulation completely without extra cost. 

We finally observe that Fig. 7(h) seems to contradict our claim that the 
combination of the two estimators always reduces the error of the better 
estimator. For the double well and the staging algorithm with L = 16 and 
32, the error of the combined estimator indeed is larger than the error for the 
kinetic estimator. The reason is that here the data for the kinetic estimator 
were in fact obtained from the simulation optimized for the virial estimator. 
As pointed out above, the choice for the optimal staging length depends on 
the estimator. The values for the two estimators and U v reported in 
Tables 1 and 2 were obtained for the respective optimal staging length j opt 
for each estimator (cp. Table 5). Clearly, the optimal staging length for 
the combined estimator would be different both from the j opt for either the 
kinetic or the virial estimator. In fact, we are dealing with an optimization 
problem of two dimensions in the parameter space of j opt and a. However, in 
practical applications the combination of the estimators would be done after 



31 



the simulation. Unfortunately, one has to decide beforehand which j opt best 
be chosen and regardless whether one might use j op t,k or j op t,v, in either case 
one would use a non-optimal j opt for one of the two estimators. In fact, if 
we would have measured the energy using the kinetic estimator using the 
optimal staging length for the virial estimator, i.e. for L = 16 with j = 6 and 
for L = 32 with j = 10 we would have obtained values of —0.9564(11) and 
of —0.9184(20) respectively. Clearly, these values are no longer superior to 
the ones obtained using the combined estimator which were also computed 
for the optimal j for the virial estimator. 

Thus, the apparent superiority of the kinetic estimator over the combined 
estimator in Fig. 7 (h) is in fact an artifact of our over-careful data analysis. 
In practice, when the combination is done after the Monte Carlo run with 
data for the same set of parameters (non-optimal for at least one of the 
estimators), the combined estimator is guaranteed to yield the best energy 
estimates with the smallest error bars. 

6 Conclusions 

Our concern in this paper was to show how energy estimation in PIMC 
simulations can be optimized by taking into account both variances and 
autocorrelation times of two standard energy estimators: the "kinetic" and 
the "virial" estimator, and by investigating their respective interplay with 
different update algorithms: the standard local Metropolis algorithm, the 
non-local multigrid V- and W-cycles, and the non-local staging algorithm. 
Let us briefly summarize the main points: 

(i) While the variance of the virial estimator depends only very weakly on 
the discretization scale, the variance of the kinetic estimator diverges 
asymptotically according to a\ = L/2f3 2 . This behaviour is indepen- 
dent of the update algorithm. 

(ii) The dynamics of the update algorithms affects the autocorrelation 
times of the standard estimators. For the Metropolis algorithm these 
diverge as L 2 but this behaviour can clearly be seen only for the virial 
estimator. In fact, the values of the autocorrelation times for the ki- 
netic estimator turned out to be much smaller than those of the virial 
estimator, and we did not see the expected L 2 divergence with the 
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parameters in our simulations. Refined non-local update algorithms 
reduce (V-cycle) or eliminate (W-cycle, staging) an L divergence of the 
autocorrelation times. For the staging algorithm we observed quite a 
strong dependence of the autocorrelation times on the length of the 
staging segment. In particular, we have shown that the optimal length 
is very sensitive to the choice of the energy estimator. For both estima- 
tors it scales, however, linearly with L. In terms of acceptance rates the 
best performance was obtained at about 55% for the virial and about 
85% — 90% for the kinetic estimator, respectively. 

(iii) The interplay between the variances of estimators and the dynamics of 
the update algorithms which affects the autocorrelation times turns out 
to be quite subtle. It furthermore also depends on the discretization 
parameter L. The kinetic estimator often has a smaller variance than 
the virial estimator for small L. Looking at the errors of the energy esti- 
mates we hence observe a crossover at which the virial estimator starts 
to win over the kinetic estimator since its variance does not increase 
with L. For the Metropolis algorithm, due to its small autocorrelation 
times for the kinetic estimator, the crossover point, however, is shifted 
to very large values of L and was not seen with the parameters in our 
simulation. 

(iv) As a simple solution to the involved interplay of the variance and the 
autocorrelation time for the two energy estimators we have introduced a 
"combined estimator." By construction this always gives more accurate 
energy estimates than the better of the two standard estimators. The 
empirically observed gain varies but realistically one may expect a gain 
of about 10% in the error for energy estimation in PIMC simulations. 
We emphasize, that this corresponds to a 20% gain in actual simulation 
time which comes at no extra cost and is very easy to implement after 
the actual simulation. 

(v) The combination of different estimators for the same physical quan- 
tity is a very general option for Monte Carlo simulations and by no 
means restricted to the use of the standard energy estimators in PIMC 
simulations. 



33 



Acknowledgments 

W.J. thanks the Deutsche Forschungsgemeinschaft for a Heisenberg fellow- 
ship. 



34 



References 

1 H. De Raedt and A. Lagendijk, Phys. Rep. 128 (1985) 233; B.J. Berne and 
D. Thirumalai, Ann. Rev. Phys. Chem. 37 (1986) 401; J.D. Doll, D.L. 
Freeman, and T.L. Beck, Adv. Chem. Phys. 78 (1990) 61. 

2 J.A. Barker, J. Chem. Phys. 70 (1979) 2914. 

3 M. Creutz and B. Freedman, Ann. Phys. 132 (1981) 427. 

4 M.F. Herman, E.J. Bruskin, and B.J. Berne, J. Chem. Phys. 76 (1982) 
5150. 

5 M. Parrinello and A. Rahman, J. Chem. Phys. 80 (1984) 860. 

6 A. Giansanti and G. Jacucci, J. Chem. Phys. 89 (1988) 7454. 

7 J.S. Cao and B.J. Berne, J. Chem. Phys. 91 (1989) 6359. 

8 W. Janke and T. Sauer, in: Proceedings of the international conference 
Path Integrals from meV to MeV, Tutzing, 1992, eds. H. Grabert, A. In- 
omata, L. Schulman, and U. Weiss (World Scientific, Singapore, 1993), 
p. 17. 

9 W. Janke and T. Sauer, Chem. Phys. Lett. 201 (1993) 499. 

10 W. Janke and T. Sauer, Chem. Phys. Lett. 263 (1996) 488. 

11 All operators will be denoted by a "hat", so that it is easy to identify 
whether a quantum statistical average (...) is to be understood in the 
operator or path-integral sense. 

12 H. Kleinert, Path Integrals in Quantum Mechanics, Statistics and Polymer 
Physics, 2nd edition (World Scientific, Singapore, 1995). 

13 Note that even though (j J2k=i f( x k)) = (f(%k)) is trivially true for any 
function f(x) due to the periodic boundary conditions in the path integral 
(5), the proper definition of estimators should always contain the summa- 
tion over the path. This becomes crucial when considering the variance 

of the estimator which involves a term ((j- J2k=i f( x k)) ) 7^ {f{%k) 2 )- 



35 



14 See, e.g., E. Merzbacher, Quantum Mechanics, 2nd edition (John Wiley, 
New York, 1970), p. 168. 

15 We would still have a valid energy estimator, though, as for any choice of 
a. 

16 With regard to the "best" estimate, given in Table 1. 

17 A.D. Sokal, Bosonic Algorithms, in Quantum Fields on the Computer, ed. 
M. Creutz (World Scientific, Singapore, 1992), p. 211. 

18 N. Madras and A.D. Sokal, J. Stat. Phys. 50 (1988) 109. 

19 W. Janke, Monte Carlo Simulations of Spin Systems, in Computational 
Physics: Selected Methods - Simple Exercises - Serious Applications, 
eds. K.H. Hoffmann and M. Schreiber (Springer, Berlin, 1996), p. 10-43. 

20 W. Janke and T. Sauer, J. Stat. Phys. 78 (1995) 759. 

21 R.G. Miller, Biometrika 61 (1974) 1; B. Efron, The Jackkmfe, the Bootstrap 
and other Resampling Plans (SIAM, Philadelphia, PA, 1982). 

22 N. Metropolis, S. Ulam, J. Am. Stat. Ass. 44 (1949) 335; N. Metropolis, A. 
W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. 
Phys. 21 (1953) 1087. 

23 U. Wolff, Phys. Lett. B288 (1992) 166. 

24 T. Sauer, Ph.D. thesis, FU Berlin (1994). 

25 J. Goodman and A.D. Sokal, Phys. Rev. Lett. 56 (1986) 1015; Phys. Rev. 
D40 (1989) 2035; G. Mack, in Nonperturbative Quantum Field Theory, 
Cargese 1987, eds. G.'t Hooft et al. (Plenum, New York, 1988), p. 309; 
G. Mack and S. Meyer, Nucl. Phys. B (Proc. Suppl.) 17 (1990) 293; D. 
Kandel, E. Domany, D. Ron, A. Brandt, and E. Loh, Jr., Phys. Rev. Lett. 
60 (1988) 1591; D. Kandel, E. Domany, and A. Brandt, Phys. Rev. B40 
(1989) 330. 

26 W. Hackbusch, Multi-Grid Methods and Applications (Springer, Berlin, 
1985). 



36 



W. Janke and T. Sauer, Phys. Lett. A197 (1995) 335. 



W. Janke, Recent Developments in Monte- Carlo Simulations of First- 
Order Phase Transitions, in Computer Simulations in Condensed Matter 
Physics VII, eds. D.P. Landau, K.K. Mon, and H.-B. Schiittler (Springer, 
Berlin, 1994), p. 29. 

W. Janke and T. Sauer, Phys. Rev. E49 (1994) 3475. 

E.L. Pollock and D.M. Ceperley, Phys. Rev. B30 (1984) 2555; M. Sprik, 
M.L. Klein, and D. Chandler, Phys. Rev. B31 (1985) 4234; ibid. B32 
(1985) 545; M.E. Tuckerman, B.J. Berne, G.J. Martyna, and M.L. Klein, 
J. Chem. Phys. 99 (1993) 2796. 

M. Takahashi and M. Imada, J. Phys. Soc. Jpn. 53 (1984) 963, 3765; X.-P. 
Li and J.Q. Broughton, J. Chem. Phys. 86 (1987) 5094. 

Recall that for the Metropolis algorithm measurements were performed 
only every n e sweep. Here we have rescaled the measured autocorrelation 
times to units of sweeps. 

For a comparison of Figs. 7(a) and (b) with Tables 1 and 2 one has to keep 
in mind that the errors given in the tables refer to the full simulation with 
n e x N m sweeps, while in the figures we have rescaled them to hypothetical 
runs with N m = 100 000 sweeps (=measurements), i.e., the numbers in 
the tables have to multiplied with ^fn~ e to obtain the data points in the 
figures. 

This is, in fact, what we expect on general theoretical grounds. 

For the Metropolis algorithm the autocorrelation times of the two estima- 
tors differ so much for L > 128 that it makes no sense to compute the 
combined estimator. 



37 







Convex Potential (CP): 


V(x) = ±x 2 


fx 4 
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8 


0.48883(57) 


0.12 


0.4891(21) 


0.43 


0.913(58) 


1.323 


0.48874(29) 


0.060 


16 


0.65741(86) 


0.13 


0.6553(23) 


0.36 


0.898(63) 


1.120 


0.65766(83) 


0.13 


32 


0.7534(18) 


0.24 


0.7552(31) 


0.42 


0.86(11) 


0.793 


0.7537(17) 


0.23 


64 


0.7864(37) 


0.47 


0.7783(58) 
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0.78(23) 


0.825 


0.7850(36) 


0.46 
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0.67(42) 
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0.8001(66) 
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0.791(11) 
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0.54(42) 


0.933 


0.8063(69) 


0.86 
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0.8169(73) 


0.89 


0.817(10) 


1.3 


0.42(39) 


0.819 


0.8169(72) 


0.88 
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0.7990(96) 


1.20 


0.794(11) 


1.4 


0.30(38) 


0.611 


0.7973(86) 


1.08 


V-cycle 
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0.48834(40) 
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0.4871(16) 
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1.253 


0.48866(25) 
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16 


0.65847(57) 
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0.6581(16) 


0.24 


0.898(62) 


0.969 


0.65846(57) 


0.087 


32 


0.7530(11) 
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0.7534(17) 
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0.857(73) 
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0.75313(81) 
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64 


0.7916(23) 


0.29 


0.7864(18) 


0.23 


0.78(14) 


0.391 


0.7885(14) 


0.18 


128 


0.8011(36) 


0.45 


0.7982(21) 


0.27 


0.67(23) 


0.261 


0.7990(18) 


0.23 


256 


0.8085(58) 


0.72 


0.8051(26) 


0.33 


0.54(35) 


0.170 


0.8057(24) 


0.30 


512 


0.8174(85) 


1.04 


0.8051(39) 


0.48 


0.42(38) 


0.192 


0.8075(34) 


0.42 
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0.807(13) 


1.61 


0.8008(56) 


0.70 


0.31(45) 


0.164 


0.8018(52) 


0.65 


W-cycle 
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0.48844(32) 


0.066 


0.4882(13) 
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0.912(42) 


1.228 


0.48850(22) 
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0.65822(53) 
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0.6589(13) 
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0.343 
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0.7963(14) 
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0.67(20) 
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0.66 


0.8045(14) 
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0.8050(14) 
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0.801(11) 


1.38 


0.8033(18) 


0.22 


0.31(39) 


0.025 


0.8032(18) 


0.22 


Staging 
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0.48786(52) 
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0.4861(19) 


0.39 


0.913(53) 


1.337 


0.48845(25) 


0.051 


16 


0.65907(67) 


0.10 


0.6594(20) 


0.30 


0.899(69) 


1.121 


0.65743(62) 


0.094 


32 


0.7544(13) 


0.17 


0.7556(21) 


0.28 


0.858(90) 


0.779 


0.7546(12) 


0.16 


64 


0.7902(20) 


0.25 


0.7889(23) 


0.29 


0.79(16) 


0.411 


0.7883(18) 


0.23 
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0.7958(38) 


0.48 


0.7980(22) 


0.28 


0.67(26) 
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0.24 
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0.8051(50) 


0.62 


0.8025(27) 
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0.54(35) 


0.101 


0.8035(26) 


0.32 
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0.7981(71) 


0.89 


0.8050(26) 


0.32 


0.41(48) 


0.076 


0.8035(25) 


0.31 
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0.799(11) 


1.4 


0.7968(22) 


0.28 


0.31(52) 


0.047 


0.7969(21) 


0.26 



Table 1: Convex Potential (CP): Measured energies using the kinetic estima- 
tor £4, the virial estimator U v , and the combined estimator U c . Also listed 
are the relative jacknife errors in percent, the cross-correlation coefficient p, 
and the parameter a opt of the optimal combined estimator. 
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-1.04822(45) 
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-1.0484(26) 
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0.88(10) 
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-1.04822(45) 
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-0.95403(81) 
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-0.95487(74) 


0.078 


32 


-0.9179(15) 


0.17 


-0.9195(19) 


0.21 
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-1.04841(65) 


0.062 


-1.0542(37) 
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0.88(13) 


1.003 


-1.04839(65) 


0.062 


16 


-0.95318(91) 


0.095 


-0.9582(35) 


0.37 
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0.916 


-0.9566(11) 


0.11 


32 


-0.9189(15) 


0.16 


-0.9157(44) 


0.48 


0.79(20) 


0.817 


-0.9179(18) 


0.20 


64 
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0.29 
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0.28 
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-0.9013(37) 
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0.51 
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0.41 


512 
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0.87 
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0.47 


0.42(35) 


0.229 


-0.9044(36) 


0.40 


1024 


-0.907(11) 


1.2 


-0.8999(44) 


0.49 


0.31(52) 


0.112 


-0.9004(41) 


0.46 



Table 2: Double Well (DW): Measured energies using the kinetic estimator 
£4, the virial estimator U v , and the combined estimator U c . Also listed are 
the relative jacknife errors in percent, the cross-correlation coefficient p, and 
the parameter a opt of the optimal combined estimator. 
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Convex Potential (CP): V(x) = ±x 2 + x i 



L 


< 


f"int,k 




Tint.v 


< 


Tint.c 


Metropolis 


8 


0.00916(13) 


1.636(51) 


0.1049(14) 


1.825(48) 


0.003440(28) 


1.053(28) 


16 


0.01942(16) 


1.758(51) 


0.1272(18) 


2.378(56) 


0.02097(15) 


1.535(32) 


32 


0.07641(58) 


2.079(45) 


0.1345(19) 


4.44(14) 


0.05683(53) 


2.879(79) 


64 


0.2249(17) 


2.619(81) 


0.1278(30) 


11.8(1.0) 


0.1569(13) 


3.40(14) 


128 


0.5418(43) 


3.30(12) 


0.1296(59) 


47.6(9.6) 


0.4055(34) 


4.07(17) 


256 


1.1760(56) 


6.03(15) 


0.1285(48) 


165(36) 


1.0247(49) 


6.47(17) 


512 


2.469(12) 


13.91(24) 


0.1311(42) 


660(110) 


1.6593(83) 


15.97(29) 


1024 


5.015(23) 


47.55(80) 


0.1280(52) 


3160(570) 


1.8888(83) 


65.9(1.9) 


V-cycle 


8 


0.00929(10) 


0.845(15) 


0.1059(11) 


1.033(16) 


0.003658(32) 


0.6988(72) 


16 


0.01965(24) 


0.746(11) 


0.1297(20) 


0.912(14) 


0.01978(26) 


0.732(11) 


32 


0.07632(46) 


0.982(17) 


0.1335(10) 


0.909(18) 


0.05284(35) 


0.824(13) 


64 


0.2248(11) 


1.163(22) 


0.1352(13) 


1.092(22) 


0.08560(71) 


0.947(18) 


128 


0.5450(32) 


1.326(19) 


0.1350(12) 


1.583(31) 


0.11181(83) 


1.399(25) 


256 


1.1847(72) 


1.460(25) 


0.1352(15) 


2.701(82) 


0.1275(12) 


2.321(65) 


512 


2.451(13) 


1.638(30) 


0.1372(22) 


5.07(21) 


0.1785(16) 


3.062(90) 


1024 


5.014(28) 


1.758(32) 


0.1345(27) 


10.80(94) 


0.2299(24) 


5.07(25) 


W-cyclc 


8 


0.009041(81) 


0.6434(93) 


0.10449(96) 


0.851(16) 


0.003796(26) 


0.5702(36) 


16 


0.01962(12) 


0.652(11) 


0.1292(12) 


0.692(12) 


0.02198(17) 


0.5485(69) 


32 


0.07590(44) 


0.902(12) 


0.13406(93) 


0.644(11) 


0.05400(36) 


0.5626(75) 


64 


0.2265(11) 


1.033(18) 


0.1362(10) 


0.6278(66) 


0.08669(57) 


0.6098(72) 


128 


0.5454(29) 


1.103(16) 


0.13335(91) 


0.660(14) 


0.10754(73) 


0.641(12) 


256 


1.1794(73) 


1.222(19) 


0.13640(85) 


0.7119(93) 


0.12236(94) 


0.690(10) 


512 


2.460(13) 


1.177(22) 


0.1357(11) 


0.781(12) 


0.12951(97) 


0.765(12) 


1024 


5.002(28) 


1.250(25) 


0.1348(12) 


0.859(16) 


0.1314(12) 


0.845(14) 


Staging 


8 


0.00900(11) 


1.398(49) 


0.1032(12) 


1.545(32) 


0.003454(26) 


0.9084(36) 


16 


0.01936(19) 


1.396(35) 


0.1275(13) 


1.518(24) 


0.02117(14) 


0.7878(39) 


32 


0.07628(45) 


1.086(21) 


0.1347(16) 


1.850(44) 


0.05588(42) 


1.287(32) 


64 


0.2270(15) 


1.114(22) 


0.1336(14) 


1.960(44) 


0.08546(75) 


1.644(41) 


128 


0.5453(29) 


1.115(23) 


0.1357(16) 


2.012(49) 


0.1083(12) 


1.885(47) 


256 


1.1881(66) 


1.110(18) 


0.1344(18) 


2.119(67) 


0.1209(15) 


2.086(64) 


512 


2.448(13) 


1.076(16) 


0.1339(16) 


2.079(53) 


0.1281(15) 


2.019(57) 


1024 


5.054(29) 


1.097(21) 


0.1319(17) 


2.048(56) 


0.1306(16) 


1.983(54) 



Table 3: Convex Potential (CP): Variances a\ N c and the integrated autocor- 
relation times Ti nti k,v,c for the kinetic (k), virial (v) and combined (c) energy 
estimator. 
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Double Well (DW): V{x) = -~\x z + 0.04x 4 



L 


< 


f"int,k 


< 


Tint.v 




Tint.c 


Metropolis 


8 


0.01282(30) 


1.833(76) 


0.3645(40) 


2.564(48) 


0.00128(30) 


1.834(76) 


16 


0.03257(41) 


1.87(13) 


0.3683(74) 


5.88(24) 


0.00322(41) 


1.89(13) 


32 


0.1050(14) 


2.68(20) 


0.364(12) 


18.7(2.1) 


0.1096(14) 


2.59(19) 


64 


0.2598(19) 


3.206(82) 


0.352(13) 


64.4(9.4) 


0.2327(17) 


3.62(11) 


128 


0.5746(42) 


3.75(20) 


0.361(11) 


243(21) 






256 


1.217(11) 


4.99(24) 


0.3500(95) 


909(67) 






512 


2.498(19) 


6.17(33) 


0.358(13) 


4407(570) 






1024 


5.048(23) 


10.68(34) 


0.365(12) 


14752(1600) 






V-cycle 


8 


0.01247(18) 


0.839(53) 


0.3682(28) 


0.752(18) 


0.00124(19) 


0.840(50) 


16 


0.03356(27) 


0.918(17) 


0.3636(29) 


0.746(12) 


0.00330(23) 


0.842(13) 


32 


0.10508(68) 


1.146(23) 


0.3590(27) 


0.784(12) 


0.00840(57) 


0.968(15) 


64 


0.2605(16) 


1.290(27) 


0.3538(26) 


0.945(16) 


0.1551(11) 


1.021(18) 


128 


0.5830(35) 


1.438(28) 


0.3575(36) 


1.338(25) 


0.2250(17) 


1.396(27) 


256 


1.2168(65) 


1.598(37) 


0.3550(41) 


2.297(61) 


0.2749(28) 


2.057(57) 


512 


2.519(14) 


1.814(41) 


0.3561(58) 


3.92(15) 


0.4293(39) 


2.498(65) 


1024 


5.037(30) 


1.790(38) 


0.3657(74) 


7.19(43) 


0.5731(65) 


3.60(14) 


W-cyclc 


8 


0.01215(13) 


0.597(20) 


0.3631(24) 


0.5478(76) 


0.00120(13) 


0.598(20) 


16 


0.03415(24) 


0.836(26) 


0.3647(25) 


0.5121(79) 


0.00340(22) 


0.709(22) 


32 


0.10389(63) 


0.987(20) 


0.3587(23) 


0.5128(74) 


0.00930(61) 


0.639(13) 


64 


0.2618(14) 


1.110(27) 


0.3573(25) 


0.5291(75) 


0.1837(13) 


0.5762(91) 


128 


0.5856(33) 


1.191(22) 


0.3630(28) 


0.5377(72) 


0.2323(17) 


0.6266(88) 


256 


1.2157(67) 


1.274(22) 


0.3566(26) 


0.5479(72) 


0.2864(20) 


0.5841(74) 


512 


2.504(14) 


1.236(24) 


0.3627(26) 


0.5650(95) 


0.3246(24) 


0.5842(92) 


1024 


5.070(28) 


1.287(23) 


0.3564(24) 


0.5920(76) 


0.3345(22) 


0.6360(84) 


Staging 


8 


0.01242(30) 


1.323(29) 


0.3631(38) 


1.779(36) 


0.01239(31) 


1.325(29) 


16 


0.03330(34) 


1.240(48) 


0.3582(41) 


2.059(44) 


0.03242(43) 


1.579(61) 


32 


0.10403(76) 


1.120(23) 


0.3628(49) 


2.182(63) 


0.08351(77) 


1.575(43) 


64 


0.2606(16) 


1.117(22) 


0.3539(48) 


2.335(64) 


0.1568(14) 


1.648(27) 


128 


0.5847(35) 


1.145(23) 


0.3582(49) 


2.406(76) 


0.2337(25) 


2.017(49) 


256 


1.2242(67) 


1.139(24) 


0.3597(52) 


2.456(87) 


0.2980(32) 


2.085(71) 


512 


2.520(14) 


1.164(20) 


0.3592(52) 


2.366(64) 


0.3473(36) 


2.095(49) 


1024 


5.067(27) 


1.150(20) 


0.3610(49) 


2.406(69) 


0.3478(42) 


2.224(59) 



Table 4: Double Well (DW): Variances <J^ VC and the integrated autocorre- 
lation times r inti k,v,c for the kinetic (k), virial (v) and combined (c) energy 
estimator. 



41 



L 


Jopt,k 


% 


Jopt,v 


% 


Convex Potential (CP) 


8 


2 


64 


2 


64 


16 


4 


49 


2 


82 


32 


4 


73 


4 


72 


64 


4 


90 


10 


68 


128 


10 


86 


24 


55 


256 


16 


90 


44 


55 


512 


24 


93 


88 


55 


1024 


56 


90 


176 


56 


Double- Well (DW) 


8 


2 


71 


2 


71 


16 


4 


63 


6 


47 


32 


4 


80 


10 


49 


64 


8 


84 


20 


52 


128 


16 


80 


40 


54 


256 


32 


82 


80 


54 


512 


64 


84 


160 


54 


1024 


128 


84 


320 


54 



Table 5: Optimal staging lengths j op t,k, jo P t, v an d acceptance rates for the 
kinetic (k) and the virial (v) estimators. 



42 







/ ■ 




c kinetic 












a combined 










§ 


s 

i i 

i 


CP, Metropolis 



100 
L 











c kinetic 






□ virial 






a combined 













m @ 

@ A 


8 * 


® 


CP, V-cycle 



100 
L 



kinetic 
□ virial 
combined 




CP, W-cycle 



100 
L 



o kinetic 

□ virial 

a combined 



CP, Staging 



100 
L 




DW, Metropolis 



100 
L 



kinetic 
□ virial 
a combined 



DW, V-cycle 



100 
L 



kinetic 
□ virial 
a combined 



100 
L 



DW, W-cycle 
1000 



o kinetic 

□ virial 

a combined 



-<s © © 



DW, Staging 



100 
L 



Figure 6: Integrated autocorrelation times r int on a logarithmic scale for 
the three energy estimators using different update algorithms for the convex 
potential (CP) and the double well (DW). Straight lines show fits of the form 
ri nt = aL z . 43 
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Figure 7: The errors A = AU computed from the variance and the autocor- 
relation times for the three estimators using different update algorithms for 
the convex potential (CP) and the double well (DW). 
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Figure Captions 



Fig. 1: Variance of the individual energy measurements using the "kinetic" 
and the "virial" estimators for the two potentials (6) and (7) at (3 = 
10. While the variance of the virial estimator is roughly constant in 
the continuum limit L — > oo, the variance of the kinetic estimator 
asymptotically diverges as = L/2/3 2 . 

Fig. 2: The optimal interpolation parameter a opt as given in eq. (26) as a 
function of the ratio k = A v /A k and the correlation coefficient p = 
AL/(A k A v ). 

Fig. 3: The reduction factor R as given in eq. (27) as a function of the ratio 
k = A v /A k and the correlation coefficient p = A kv /(A k A v ). 

Fig. 4: Relative statistical error AU C /U C of the linear combination U c of 
the two energy estimators for the convex potential as a function of 
the interpolation parameter a (using the W-cycle update algorithm). 
Data symbols denote jackknife averages over 100 blocks, and the solid 
lines are computed according to the theoretical prediction (23), using 
as input only the (jackknife) variances of the virial and the kinetic 
estimator (corresponding to a — 0, resp. a — 1) and the (jackknife) 
covariance of the two estimators. The arrows indicate the values of 
optimal a according to eq. (24). 

Fig. 5: (a) The autocorrelation function A(j) on a logarithmic scale and 
(b) the integrated autocorrelation time T int {k) of the virial estimator as 
obtained with the V-cycle multigrid algorithm for the convex potential 
(CP) at (5 — 10 and L = 512. Solid lines show fits according to 
eq. (31) resp. (33). The asymptotic value of r in t(A;) quoted in Table 3 
is r intiV = 5.07(21). 

Fig. 6: Integrated autocorrelation times r int on a logarithmic scale for the 
three energy estimators using different update algorithms for the convex 
potential (CP) and the double well (DW). Straight lines show fits of 
the form r int = aL z . 
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7: The errors A = AU computed from the variance and the autocorre- 
lation times for the three estimators using different update algorithms 
for the convex potential (CP) and the double well (DW). 
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